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We propose a novel coarse graining tensor renormalization group method based on the higher- 
order singular value decomposition. This method provides an accurate but low computational cost 
technique for studying both classical and quantum lattice models in two- or three-dimensions. We 
have demonstrated this method using the Ising model on the square and cubic lattices. By keeping 
up to 16 bond basis states, we obtain by far the most accurate numerical renormalization group 
results for the 3D Ising model. We have also applied the method to study the ground state as well 
as finite temperature properties for the two-dimensional quantum transverse Ising model and obtain 
the results which are consistent with published data. 
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I. INTRODUCTION 

The simulation of two or higher dimensional quan- 
tum lattice models remains a great challenge. This has 
stimulated great interest on the investigation of renor- 
malization group (RG) methods for the tensor-network 
states-^—. The use of the tensor-network state as a vari- 
ational wave function for the classical lattice model was 
first considered by Nishino and coworkers^—. They, 
and recently Garcia-Saez et. a/ 17 , proposed a number 
of RG approaches to study the thermodynamic proper- 
ties of the Ising and other models. However, due to the 
heavy computational cost, the maximal truncated tensor 
dimension, D, that can be handled with their methods is 
small (between 2 and 5) in 3D, and consequently the ac- 
curacy of the results they obtained is low in comparison 
with the Monte Carlo ones. 

In 2007, Levin and Nave^ proposed a coarse grained 
tensor renormalization group (TRG) method for study- 
ing two dimensional (2D) classical models based on the 
singular value decompostion (SVD) of matrix. Later 
we proposed a second renormalization group (SRG) 
method 7,8 to globally optimize the truncation scheme 
and improve significantly the accuracy of the TRG. The 
application of these methods in classical and quantum 
lattice models has achieved great success^£i£~— . How- 
ever, it is difficult to extend these methods to 3D, not 
just due to the increase of the order of local tensors, but 
also due to the change of lattice topology in the coarse 
graining process^. 

In this paper, we introduce a novel coarse graining 
TRG method based on the higher-order singular value 
decomposition (HOSVD )2£ to study physical properties 
of 2D or 3D lattice models. We will first discuss about a 
simple TRG method based on the HOSVD (abbreviated 
as HOTRG hereafter), and then discuss about a more so- 
phisticated method that incorporates the second renor- 



malization effect of environment tensors to the HOTRG. 
The HOSVD takes into account more accurately the in- 
terplay between different components of a tensor. It pro- 
vides a better scheme to truncate a local tensor than the 
SVD. 

This paper is arranged as follows. In Sec. [Hi a de- 
tailed introduction to the HOTRG and HOSRG for the 
2D statistical lattice models is given. We have taken 
the 2D Ising model to show how accurate the HOTRG 
and HOSRG can be in comparison with other methods. 
In Sec. HB we have extended the HOTRG and HOSRG 
to the 3D statistical lattice models. For the 3D Ising 
model, we have obtained by far the most accurate nu- 
merical renormalization group results for the 3D Ising 
model. Our accuracy is comparable with the best Monte 
Carlo results. In Sec.QVl we have applied HOTRG to 2D 
quantum lattice models. Our preliminary results show 
that the HOTRG provides a powerful tool for studying 
the ground state and thermodynamic properties of 2D 
quantum lattice models. A summary is given in Sec. |V] 



II. TWO-DIMENSIONAL SYSTEMS 
A. HOTRG 

Let us start by taking the Ising model, 

h = -J24< (1) 

m 

as an example to show how the method works in 2D 
first. a\ is the pauli matrix at site i. An extension of 
the method to 3D will be described later. The partition 
function of the 2D Ising model can be represented as a 
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FIG. 1: (a) A HOTRG contraction of the tensor network 
state along the y axis on the square lattice, (b) Steps of 
contraction and renormalization of two local tensors. The 
initial tensor T (0) = T. 



translation invariant tensor network stated, 



(2) 



where i runs over all the lattice sites and Tr is to sum 
over all bond indices, and the local tensor T is defined at 
each lattice site as shown in Fig.QJa), 



a 

where W is a 2 x 2 matrix defined by 



(3) 



/ y /cosh(l/T) , y/ sinh(l/T) \ 

^ Vcosh(l/T), -Vsinh(l/T) J ' lJ 

and T is the temperature. 

To coarse grain, we contract the lattice alternatively 
along the horizontal (x-axis) and vertical (y-axis) direc- 
tions. This scheme of coarse graining is simple to imple- 
ment. Fig. [IJa), as an example, shows how the contrac- 
tion along the y-axis is done. At each step, two sites are 
contracted into a single site in the coarse grained lattice 
(Fig. Hlb)), an d the lattice size is reduced by a factor of 
2. 

The contracted tensor at each coarse grained lattice 
site is defined by 



M 



Ey,(n) y,(n) 



(5) 



where x = x± ® X2, x' = x'^ ® x' 2 , and the superscript n 
denotes the n'th iteration. The bond dimension of 
along the x-axis is the square of the corresponding bond 



dimension of T^ n \ To truncate into a lower rank 

tensor, we first do a HOSVD for this tensor— 
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(6) 



where U's are the unitary matrices. S is the core tensor 
of M("), which possesses the following properties for any 
index, say index j: 
(1) all orthogonality, 



(S :J ,,\S : , r ..) = 0, itj^f, 



(7) 



where (S : j i:i : | S : j' i: , : ) is the inner-product of these two 
sub-tensors. 

(2) pseudo-diagonal, 



\S : , j .,\>\S:,r-,:\, Xj<f, 



where IS- 



is the norm of this sub-tensor which is the 



square root of all elements' square sum. These norms 
play a similar role as the singular values of a matrix. 

In M^ n \ the two vertical bonds, y and y', do not need 
to be renormalizcd. Thus in the practical calculation, 
U u and U D are not needed to be determined. Moreover, 
the right bond of is linked directly to the left bond 
of an identical tensor on the right neighboring site, thus 
to truncate any one of the horizontal bonds of M W will 
automatically truncate the other horizontal bond. The 
truncation can be done by comparing the values of 
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(8) 



(9) 



If £1 < £2, we truncate the first dimension of S or the 
second dimension of U L to D. Otherwise, we truncate 
the second dimension of S or the second dimension of 
U R to D. This kind of truncation scheme provides a 
simple and optimal approximation to minimize the trun- 
cation erro r 21 ' 22 . It has been successfully applied to many 
fields such as data compression, image processing, pat- 
tern recognition, and eto2£. 

After the truncation, we can update the local tensor 
using the following formula 



T 



(n+l) 



xx'yy' / _j 
ij 



ijyy jx' 



(10) 



where = U L (or U R ) if E\ is smaller (or larger) 

than £2- 

The above HOTRG calculation can be repeated iter- 
atively until the free energy and other physical quanti- 
ties calculated are converged. The cost of the calculation 
scales as D 7 in the computer time and D 4 in the memory 
space. This is comparable with the cost of TRG££. 
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FIG. 2: (color online) Graphical representation of Eq. (|13[) 
for determining the environment tensor Ej^j ^ from E^™^ 
in the backward iteration. 

The key step in the above HOTRG iteration is to deter- 
mine the four unitary matrices on the right hand side of 
Eq. ([6]). In our calculation, we determine these matrices 
by taking the singular value decomposition of matrices. 
As an example, let us consider how to evaluate U L . We 
first convert M xx i yy i into a matrix M' x x , , 

lvl x,x'yy' ~ lvl xx'yy' 

with the first index x as the row index and the rest indices 
(a/, y, y') as the column index of this matrix. Then from 
the theory of HOSVD, we know that U L is equal to the 
left unitary matrix of M under the singular value decom- 
position. Thus U L can be simply determined from the 
canonical transformation of the unitary matrix M' 'M'* 

M'M ,jf = U L K L {U L )\ (11) 

where A L is the eigenvalue of M'M'^. Furthermore, it 
can be shown that 

|S i): , :>: | 2 = Af. (12) 
The cost for evaluating these U- matrices scales with D 6 . 

B. HOSRG 

The HOTRG is a local optimization method. It mini- 
mizes the error in the truncation of a local tensor. How- 
ever, it ignores the renormalization effect of environment. 
To develop a global optimization method, it is necessary 
to consider the environment contribution in the renor- 
malization of local tensors. In Refsji^, we proposed a 
SRG appraoch to incorporate the environment contri- 
bution in the optimization of local tensors. This kind 
of SRG approach can be also used to improve the per- 
formance of HOTRG, which leads to a global optimized 
HOTRG method, referred as HOSRG below. 

The HOSRG follows the same coarse graining steps 
as in the HOTRG. However, at each step, one needs to 



calculate a bond density matrix defined on a bond whose 
basis space will be truncated. This bond density matrix 
is defined by tracing out all environment tensors. 

The SRG introduced in Ref<££ is an infinite lattice al- 
gorithm. In the calculation of the bond density matrix, 
the size of the environment is always assumed to be infi- 
nite. At each step of coarse graining, a combined forward 
and backward iteration is performed. In the forward iter- 
ation, the TRG is applied to determine all transformation 
matrices and local tensors. A backward iteration is then 
performed to determine the bond density matrix. This 
scheme can be readily extended to the HOSRG. We have 
done this kind of calculation and find that it does provide 
much more accurate results than the HOTRG. 

Similar as in the DMRG, one can also introduce a fi- 
nite lattice algorithm to perform the HOSRG calcula- 
tion. In this case, the whole size of the system is fixed 
and the number of tensors in the environment is re- 
duced at each step of coarse graining. Again the second 
renormalization effect of the environment is handled by 
performing forward-backward iterations. But now this 
kind of forward-backward iterations can be repeated for 
many times, similar as to do a finite size sweeping in the 
DMRG. This provides a self-consistent approach to treat 
the system as well as environment tensors. It can further 
improve the accuracy of the HOSRG. Below we give an 
introduction to this method. 

In the first round of forward iteration, we carry out a 
standard HOTRG calculation to determine iteratively all 
the transformation matrices U^ n > and local tensors T", 
This iteration ends when the system reaches a desired 
size, say 2 N lattice points with ./V = 30 ~ 50. One then 
carries out a backward iteration to calculate the environ- 
ment tensor £K n J iteratively, starting from E^ N+1 ' which 
is set to be an unit tensor. The iteration formula for 
determining j s given by 

p(") _ p(n+lWn) T j(n+l)jj(n+l) „s 

■^kajtii ~ *jkl 1 i 2j2 al u i 1 i 2 ,i U hh,j ■ \ °> 

A graphical representation of this equation is shown in 
Fig. [2j This backward iteration is terminated after E^ 
is determined. 

The above forward-backward iteration determines all 
coarse grained system and environment tensors that are 
needed for carrying out the second renormalization cal- 
culation. From these tensors we can perform another for- 
ward and backward iteration to improve their accuracy. 
But starting from this round of iteration, the HOTRG is 
no longer needed. Instead, at each step of coarse grain- 
ing we evaluate the bond density matrix p zw ,xy using the 
following formula 

(n) p (n+2) Tr (n+l) TT (n+l) TT (n+2) TT (n+2) 
Pzw.xy - a ijkl u ni 2 i U hhj u k l k 2 k u l 1 l 2 l 

rp(n) rp(n) rp{n) rp(n) (- 14 % 
i\xk\a i 2 yal± zj±k 2 b wj 2 bl 2 ' ^ ' 

The repeated index summation is assumed. A graphical 
representation of this formula is shown in Fig. [3] 
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FIG. 3: (color online) Graphical representation of Eq. (|14[l for 
determining the bond density matrix pi$ tX y through E 1 :™^. 
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FIG. 4: (color online) Comparison of the relative errors of 
free energy with respect to the exact results for the 2D Ising 
model obtained by various methods with D = 24. The critical 
temperature T c = 2/ln(l + \/2). 
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FIG. 5: (a) A HOTRG coarse graining step along the z-axis on 
the cubic lattice, (b) Steps of contraction and renormalization 
of two local tensors. 



result is already less than 10 7 even at the critical tem- 
perature, much more accurate than the TRG result^. 
The HOSRG also performs better than the SRG. But the 
difference in the results obtained by these two methods is 
relatively small around the critical point. The HOTRG is 
less accurate than the two SRG methods, but it is compu- 
tationally economic. The difference between TRG/SRG 
and HOTRG/HOSRG lies mainly in the basis truncation 
scheme. The former is based on the SVD, while the latter 
is based on the HOSVD. The above comparison indicates 
that the HOSVD scheme works better. 



To diagonalize this bond density matrbssi 

p (n) _ J7(n+1) A (u^ n+1 A , (15) 

we can find its eigenpair, (A, [/("+ 1 )). Same as in the den- 
sity matrix renormalization group^i, the eigenvalues of 
this density matrix determine the probabilities of the cor- 
responding eigenvectors in the virtual bond basis space. 
By keeping the largest D eigenvalues of A and the cor- 
responding eigenvectors of U^ n+1 \ one can update the 
local tensor using Eq. (p~0|) . After finishing this 

forward iteration, we can take a backward iteration to up- 
date all environment tensors with Eq. (|13p . This forward- 
backward iteration is then repeated until all system and 
environment tensors are converged. 

Fig. [4] compares the relative errors of free energy with 
respect to the rigorous solution^ for the 2D Ising model 
obtained with four different methods. By keeping just 
24 states, we find that the relative error of the HOTRG 



III. THREE-DIMENSIONAL SYSTEMS 

The above HOTRG and HOSRG methods can be read- 
ily extended to 3D. This is an advantage of the coarse 
graining scheme proposed here. On the cubic lattice, a 
full cycle of lattice contraction needs to be done in three 
steps, along the x-axis, y-axis, and z-axis, respectively. 
At each step, two neighboring tensors will be combined 
to form a single coarse grained tensor and the lattice size 
is reduced by a factor of 2. 

As an example, Fig. [S] shows how the tensors are con- 
tracted along the z-axis. The HOSVD of the coarse 
grained local tensor (Fig. [H[b)) can be similarly done as 
for the 2D case. But the local tensor now has six bond 
indices and a HOSVD for a higher order tensor should be 
done. Moreover, the basis spaces for both the x-axis and 
y-axis bonds need to be renormalized. Thus we should 
determine from the core tensor and the unitary matrices 
of M^ n > not only the transformation matrix for the x- 
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FIG. 6: (Color online) Graphical representation for the de- 



termination of the environment tensor E. 
in 3D. 
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direction bonds, EA n ), but also the transformation matrix 
for the y-direction bonds, V^ n \ After that the dimen- 
sions for both x-axis and y-axis bonds are truncated and 
the local tensor is updated using and^W. The con- 
traction and renormalization of tensors along other two 
directions can be similarly done. This three-step itera- 
tion can then be repeated until the results are converged. 

After the above HOTRG iteration, one can also do a 
backward iteration to evaluate the environment tensors 
and carry out the HOSRG calculation in 3D. A graph- 
ical representation for iteratively determining the envi- 
ronment tensor in this backward iteration is shown in 
Fig. [51 A series of forward-backward iterations is then 
performed to take into account the second renormaliza- 
tion effect of the environment to the coarse grained ten- 
sors. In the subsequent forward iterations, we evaluate 
and diagonalize the bond density matrix (see Fig. [7]) and 
update the coarse grained tensors. The environment ten- 
sors are evaluated again in the backward iteration. 

In the 3D calculation, the computational time scales 
with D 11 and the memory scales with D 6 . This cost in 
the computational resource is significantly smaller than 
in other 3D numerical RG methods ^ 17 ' 19 . We have 
studied the 3D Ising model using the HOTRG for D up 
to 16. 

The temperature dependence of the internal energy U 
and the specific heat C for the 3D Ising model obtained 
by the HOTRG with D = 14 is shown in Fig. [5] and 
compared with the Monte Carlo result 27 . Our result of 
the specific heat agrees with the Monte Carlo one. At the 
critical temperature, T c = 4.511544, the internal energy 
is found to be U c = -0.995592 for D = 14. This value 
of U C) as shown in Table U also agrees well with other 
published data. 

/From the temperature dependence of the specific heat 
around the critical point, one can estimate the critical 
exponent of the specific heat with the formula 



where t = |1 — T/T c \. However, as the specific heat data 
are obtained simply from the numerical derivative of the 
internal energy, the accuracy of the specific heat data 
is much less than that of the internal energy, especially 
around the critical point. This causes a big error in the 
determination of the exponent a with the above formula. 
This problem can be solved by directly evaluating this ex- 
ponent from the temperature dependence of the internal 
energy. From the temperature integration of the specific 
heat, it is simple to show that the internal energy should 
exhibit the following critical behavior 



U = U c + at + bt 



l-a 



(17) 



where a and b are unknown parameters which can be 
determined by fitting. 

Fig. [§] shows the fitting curves for the internal energy 
around the critical point obtained with Eq. (fT7|). The 
critical exponent is found to be a = 0.1023 and 0.1137 
for the temperature higher and lower than the critical 
value, respectively. These values of the critical exponent 
are consistent with the result obtained from the series 
expansion 2 ^, 0.104, and the Monte Carlo calculation 29 , 
0.111. 

Fig. [TU] shows the temperature dependence of the spon- 
taneous magnetization M obtained by the HOTRG with 
D=14. Our data agree well with the Monte Carlo 
results 35 . From the singular behavior of M, we find that 
the critical temperature T c — 4.511615 for D = 14. Fur- 
thermore, by fitting the data of M in the critical regime 
with the formula 



M ~ f, 



(18) 



we find that the exponent 7 = 0.3295, consistent with the 
Monte Carlo 2 ^ (0.3262) and series expansion^ (0.3265) 
results. 
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FIG. 7: (Color online) Graphical representation for the de- 
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FIG. 8: (Color online) The internal energy and the specific 
heat for the 3D Ising model obtained by the HOTRG with D 
= 14. The Monte Carlo result (black curve) obtained from an 
empirical fit formula given in Ref. [53] is shown for compari- 
son. 



Fig. [TT] shows the critical temperature T c determined 
from the singular points of the internal energy as well 
as the magnetization for D up to 16. The values of T c 
obtained from these two quantities agree with each other. 
For D = 16, Tc obtained from the internal energy and the 
magnetization are 4.511544 and 4.511546, respectively. 
The relative difference is less than 10~ 6 . But T c docs 
not vary monotonically with D. It becomes converged 
only when D > 13, indicating the importance of keeping 
a large D in the 3D TRG calculation. The error in T c , 
estimated from the difference between the values of T c for 
D = 15 and D = 16, is also less than 10~ 6 . Our results 
agree with the Monte Carlo data^Zr— . 

The above discussion indicates that the HOTRG works 
very well in 3D. The accuracy of the results can be fur- 
ther improved by applying the HOSRG. However, the 
HOSRG calculation costs much more CPU time. A thor- 
ough study with the HOSRG on the 3D Ising model is 
still in progress and the results will be published sepa- 



TABLE I: Comparison of the internal energy at the critical 
temperature U c for the 3D Ising model obtained by different 
methods. 



method 



U c 



HOTRG (D = 16) 
Series expansion 30 
Series expansion 31 
Series expansion 32 
Monte Carlo^ 
Monte CarlrjS 
Monte Carlc^i 



-0.990842(3) 

-0.991(1) 

-0.9902(1) 

-0.99218(15) 

-0.990604(4) 

-0.9904(8) 

-0.990(4) 



FIG. 9: (Color online) The internal energy (D = 14) and its 
fitting curves with Eq. (|17[) around the critical point for the 
3D Ising model, a is the critical exponent for the specific 
heat. 



rately. 



IV. GROUND STATE AND 
THERMODYNAMICS OF 2D QUANTUM 
LATTICE MODELS 

A d-dimensional quantum lattice model is equivalent 
to a (d+ l)-dimensional classical model, the HOTRG and 
HOSRG methods above introduced can be also extended 
to study the ground state and thermodynamic proper- 
ties of ti-dimensional quantum lattice models. For one- 
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FIG. 10: (color online) Temperature dependence of the mag- 
netization for the 3D Ising model (D = 14). The Monte Carlo 
result is from Ref. [3^]. Inset: Logarithmic plot the magneti- 
zation around the critical point. The slope of the fitting curve 
gives the critical exponent of the magnetization, 7 = 0.3295. 
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FIG. 11: (Color online) The critical temperature T c as a func- 
tion of the bond dimension D for the 3D Ising model obtained 
from the internal energy (U) and magnetization (M), respec- 
tively. 



dimensional quantum lattice models, there are already 
many mature methods for studying the ground state as 
well as the thermodynamic properties. For example, the 
ground state can be studied by the DMRG - and the 
thermodynamics can be studied by the quantum trans- 
fer matrix renormalization group (TMRG ) ' . Here we 
will only discuss how to apply the HOTRG/HOSRG to 
a 2D quantum lattice model. 

As an example, we will take the 2D quantum Ising 
model with a transverse field to show how these methods 
work. The Hamiltonian of this model is defined by 

H = -J2°>l-hJ2°l- ( 19 ) 

(ij) ' 



TABLE II: Comparison of the critical point T c for the 3D 
Ising model obtained by different methods. 



method 


T c 


HOTRG(D=16, from U) 


4.511544 


HOTRG(D=16, from M) 


4.511546 


Monte Carlo^i 


4.511523 


Monte Carlo^ 


4.511525 


Monte Carlo^ 


4.511516 


Monte Carlo^ 


4.511528 


Series expansion 40 


4.511536 


CTMRG 1 ^ 


4.5788 


TPVAH 


4.5704 


CTMRG 14 


4.5393 


TPVA±£ 


4.554 


Algebraic variation 41 


4.547 



FIG. 12: (Color online) The ground state energy Eq, the mag- 
netization M x = (<Tx) and Mz = (o" z ) versus the applied field 
h for the 2D quantum Ising model obtained by the HOTRG 
with D = 14. 

We start by representing the partition function of this 
model as a tensor-network model in the 2+1 dimensions. 
By using the Trotter-Suzuki decomposition formula, we 
can express the partition function as 19 

Z = Tre- 0H « Tr [e^'e^*-] L + 0(r 2 ) (20) 
where 

H z = -5>X, (21) 

H x = -hJ24- (22) 

i 

(3 = Lt is the inverse temperature and t is a small Trotter 
parameter. This partition function can be also expressed 
as a product of evolution matrix V 

Z = TrV L , (23) 

where 



TABLE III: The critical field h c for the 2D quantum Ising 
model with transverse field obtained by different methods. 



method 


h c 


HOTRG(D=14) 


3.0439 


Monte Cark42 


3.044 


Series Expansion 4 ^ 


3.044 


iPEPS^ 


3.06 


VDMA±£ 


3.2 


terg£ 


3.08 


iPEPS^ 


3.04 


CTM 19 


3.14 
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FIG. 13: (Color online) Temperature dependence of the in- 
ternal energy in three different fields. The solid dots and 
open circles are obtained with the 3D coarse graining HOTRG 
with D — 14 and the imaginary time evolution approach with 
D = 24, respectively. 

is the evolution operator between two neighboring Trot- 
ter layers. To insert the complete basis set between any 
two of the exponential terms on the right hand side of 
Eq. ([20]) , it is straightforward to show that V can be ex- 
pressed as a product of local tensors. From this, we can 
express the partition function as a 3D tensor-network 
model 

Z a Tr IJ^WWi + °( t2 )> ( 25 ) 

i 

where the six-indexed local tensor is defined by 

T lrfbtlv = £ W' al W' !Tr W' !Tf W' (7b P^P.u, (26) 

with 

w , = Vcoshr, V sinh r j 
I Vcoshr, — Vsinhr I ' 

and 

1 / e Th / 2 , e r Th / 2 \ 
P= V2\ e7 h l\ -e-^/z ) ■ (28) 

The operator V governs the basis state evolution along 
the Trotter direction. The matrix element of V between 
two sets of basis states, {{vi}\ and is defined by 

tracing out all spatial indices of local tensors T in a given 
Trotter layer 

({viWHlH}) = Tr'l[T hrifihfliVi , (29) 

i 

where Tr' is to trace out all spatial indices only. 

The transverse field h appears only in the P-matrix. 
When h = 0, it is simple to show that PP^ — I. In 




T 

FIG. 14: (Color online) The specific heat versus temperature 
obtained by the imaginary time evolution approach with D = 
24 for the 2D quantum Ising model. 



this case, the Hamiltonian returns to the 2D classical 
Ising model and Eq. (|2"0)) is reduced to a product of L 2D 
tensor-network model for the Ising model with an inverse 
temperature r. 

At zero temperature, L oo, the partition function is 
a product of infinite many tensors along all three direc- 
tions. The 3D HOTRG or HOSRG method introduced 
previously can be directly applied to study physical prop- 
erties of the ground state. In this case, there is no need 
to evaluate the ground state wavefunction. This is an ad- 
vantage of this approach in comparison with other meth- 
ods which are based on the tensor-network representation 
of ground state wavefunction. 

Fig. [T2l shows the field dependence of the ground state 
energy Eq, the transverse magnetization M x and the 
longitudinal magnetization M z for the transverse Ising 
model obtained by the HOTRG with D = 14 and r = 
0.01. In agreement with other calculations, we find that 
this model exhibits a phase transition at a finite field. 
The critical field is found to be h c = 3.0439, consistent 
with other published data, as shown in Table Mil 

At finite temperature, the lattice dimension of the ten- 
sor network model is finite along the imaginary time di- 
rection. Now two approaches can be used to evaluate 
the partition function. The first is to follow the steps of 
3D HOTRG to contract the tensors alternatively along 
the three directions. The contraction along the Trotter 
direction is terminated if all bond variables along that 
direction are contracted. The iteration is then carried 
out purely along the two spatial directions as for a pure 
2D classical model. This is an accurate and efficient ap- 
proach for evaluating physical quantities. However, the 
number of temperature points that can be studied with 
this approach is quite limited for a given t, since the tem- 
perature is reduced by a factor of 2 at each contraction 
along the Trotter direction. 

The second approach is to do the imaginary time evo- 
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FIG. 15: (Color online) Transverse magnetization as a func- 
tion of temperature for the 2D quantum Ising model. The 
solid dots and open circles are obtained with the 3D coarse 
graining HOTRG with D = 14 and the imaginary time evo- 
lution approach with D — 24, respectively. 
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FIG. 16: (Color online) Longitudinal magnetization as a func- 
tion of temperature for the 2D quantum Ising model. The 
solid dots and open circles are obtained with the 3D coarse 
graining HOTRG with D = 14 and the imaginary time evo- 
lution approach with D = 24, respectively. 



lution. In this case, one starts from a one- Trotter layer 
tensor product system, whose tensor operator is denned 
by V . At each step of evolution, one more Trotter layer 
(i.e. one V operator) is added to the system and the trun- 
cation for the spatial bond dimension is done using the 
HOTRG. By performing this imaginary time evolution 
up to certain temperature, physical quantities are then 
evaluated at that temperature by tracing out all bond 
variables with the HOTRG. The inverse temperature in- 
creases linearly with the evolution number. This allows 
us to collect more temperature points. However, as the 
truncation error is accumulated with the evolution, the 
results such obtained would become less and less accurate 
with decreasing temperature. Thus this approach should 
be used only for evaluating thermodynamic quantities in 
high temperatures. A similar imaginary time evolution 
approach was recently proposed for evaluating thermo- 
dynamic quantities by Ran et al. 4 -- based on the bond 
entanglement mean field approach proposed in RefA 

In our calculation, both approaches have been used. 
The first approach (i.e. the 3D coarse graining HOTRG) 
is more accurate than the second one (i.e. the imagi- 
nary time evolution approach) especially in low temper- 
atures. It is used for collecting the low temperature data. 
The second approach allows more temperature points to 
be evaluated and is applied to evaluate thermodynamic 
quantities in high temperatures. 

Figs. n~5¥T51 show the internal energy, the specific heat, 
the transverse and longitudinal magnetization as a func- 
tion of temperature for the 2D quantum Ising model 
in three different applied fields, respectively. The solid 
and open circles are results obtained by the 3D HOTRG 
(D = 14) and the imaginary time evolution (D = 24) ap- 
proaches, respectively. The results obtained with these 
two approaches agree with each other in the interme- 
diated temperature range. In higher temperature, the 



truncation error of the imaginary time evolution is rela- 
tively small and the results obtained with this method is 
more accurate. In low temperature the results obtained 
by the 3D HOTRG are much more accurate than those 
obtained by the imaginary time evolution, since the er- 
ror is accumulated at every step of coarse graining or 
time evolution and the 3D HOTRG can reach low tem- 
peratures in a much fewer coarse graining steps. In all 
calculations, the Trotter step is set to r = 0.01. 

When h < h c , as the ground state is spontaneously 
symmetry broken with a finite magnetic order, a finite 
temperature phase transition is expected. Such phase 
transition is confirmed by our calculation, which can be 
clearly seen from the temperature dependence of the lon- 
gitudinal magnetization M z (Fig. \W\i . The critical tem- 
perature decreases with increasing field h for h < h c . The 
specific heat is obtained from the numerical derivative of 
the internal energy shown in Fig. 1131 It shows a singular 
behavior around the critical point for h — 2. 

The above discussion indicates that the HOTRG (in- 
cluding the imaginary time evolution approach) provides 
a simple and powerful method for studying the ground 
state and finite temperature properties of 2D quantum 
lattice models. By applying the HOSRG to the above 
calculation, we believe that the accuracy of results can 
be further improved. But this takes a longer time to do 
the calculation. 



V. SUMMARY 

In summary, we have proposed the HOTRG and 
HOSRG methods for studying physical properties of clas- 
sical or quantum lattice models in 2D or 3D. By compari- 
son with the exact solution of the 2D Ising model, we have 
shown that the simple HOTRG calculation can already 
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produce very accurate numerical results. The HOSRG 
takes into account the second renormalization effect of 
the environment tensors, and it can significantly improve 
the accuracy of the HOTRG. These method allow us to 
retain an unprecedentedly high bond dimension in the 
basis truncation and yield by far the most accurate nu- 
merical RG results for the 3D Ising model. We have also 
applied the HOTRG to study the ground state and ther- 
modynamic properties of the 2D quantum Ising model 
with a transverse magnetic field. Our results agree well 
with all published data. Symmetry or good quantum 
number of the tensor network state for these models can 
be used to reduce the computational and storage cost. 
This will allow us to retain more basis states to reduce 
the truncation error and further improve the accuracy of 
results. 



In this work, we have taken two translation invariant 
tensor network models, namely the Ising model and the 
transverse Ising model, as examples to show how the 
methods work. However, it should be emphasized that 
our methods work more generally. They can be extended 
to a system which is translation invariant by shifting two 
or more lattice sites, or even to a random system, such as 
a spin glass model. By combining with the local update 
method of quantum tensor product wavefunction intro- 
duced in Ref. 0, one can also use this method to study 
ground state properties of 3D quantum lattice models. 

We wish to thank H.W.J. Blote for sending us the 
Monte Carlo result shown in Fig. [8] We are indebted 
to Prof. T. Nishino for helpful discussion. This work 
was supported by NSFC (Nos. 10934008 and 10874215) 
and MOST 973 Project (No. 2011CB309703). 
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